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We present a novel method for combining the analog and photon-counting measurements of lidar 
transient recorders into reconstructed photon returns. The method takes into account the statistical 
properties of the two measurement modes and estimates the most likely number of arriving photons 
and the most likely values of acquisition parameters describing the two measurement modes. It 
extends and improves the standard combining ("gluing") methods and does not rely on any ad hoc 
definitions of the overlap region nor on any background subtraction methods. 
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PACS numbers: 010.3640, 030.5260, 040.5250. 



I. INTRODUCTION 

In order to extend the total dynamic range of mea- 
surements the same back-scattered return signal in mod- 
ern lidar acquisition systems (a.k.a. transient recorders) 
is usually sampled with two distinct methods: a fast 
analog-to-digital converter and a photon counting unit 
[1]. The two range-resolved traces are then combined 
("glued") by first correcting the dead-time effects of the 
photon counting mode [2] and then by calibrating (fit- 
ting) the analog trace to the photon counts in some suit- 
able photon-counting rate interval [3-5]. The final step 
in the construction of the "glued" trace involves choos- 
ing a suitable signal size above which only rescaled ana- 
log values are considered and below which only the 
photon-counting trace is used. The general usability of 
such "gluing" methods is hampered by several intrin- 
sic weaknesses: the "background" is usually subtracted 
from both measurement modes whereas it could be re- 
tained and used as additional information; a large va- 
riety of regressions and x^ minimizations are used in 
the calibration of the analog signal; arbitrary and not 
well defined photon-counting rate fitting intervals are 
imposed in order to stabilize the former minimizations; 
photon-counting nonlinearity is usually corrected only 
with manufacturer-supplied dead-time values [6]; the 
unused half of the measured data is simply discarded 
and the actual position of the crossover between the ana- 
log and the photon-counting trace is not selected by any 
strict rules. 



II. MEASUREMENT 

This new method is based on the elementary observa- 
tion that the same input signal is evaluated by two differ- 
ent measurement techniques. Employing simple models 
of these two measurement processes we will make use 



of all available data to construct a new composite trace of 
arriving photon numbers in such a way that the new hy- 
pothetical number of photons would in fact most likely 
produce the two actually measured traces. 

For the sake of clarity we will keep the measurement 
models simple enough to illustrate the main properties 
of the method. Nevertheless, the procedure is highly 
flexible and if greater levels of detail are required, then 
the descriptions in Eq. (l)-(4) can be simply updated 
with more complex descriptions of the two measure- 
ment processes. 



A. Analog signal 

In a typical transient recorder, the analog signal is con- 
structed by integrating the current from a photomulti- 
plier (PMT) in a sampling time At which is then dis- 
cretized by an analog-to-digital converter into the ana- 
log lidar trace. Since the PMT is a fairly linear sensor 
of arrived photons p, we can thus describe their trans- 
formation into an analog signal a with a simple linear 
transformation 



a = A{p) = ccp + ^ 



(1) 



where a is related to the PMT and amplifier gain, con- 
verting the number of incoming photons into the adc 
units. j6 is a small hardware-imposed offset (baseline) 
which enables the detection of a possible signal under- 
shoot and post-measurement determination of a true 
zero. Given the number of the input photons p, the vari- 
ance, V[fl], of the analog signal resulting from the noise 
in this chain of electronics can be, at least for small sig- 
nals, safely modeled as being constant. 



V[a] = 7^ 



(2) 
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and is expressed in units of adc^. For larger signals the 
analog variance in Eq. (2) acquires an additional signal- 
size dependent Poisson term which we can neglect for 
reasons given later. 



B. Photon-counting mode 

In t5rpical photon-counting modules of modern tran- 
sient recorders, the input photons p are recorded by 
counters with predominantly non-extending dead-time 
T. These types of counters are also referred to as cumu- 
lative or non-paralyzable counters. 

For such counters^ the mean number of counts mina 
sampling time At can be expressed as 



m 



C{p) 



1 + Sp' 



(3) 



where S is the fraction of dead-time vs. sampling time, 
S = t/ At. The variance of the photon count is 



V[m] = V,{p), 



(4) 



where V^ is a nontrivial function for the variance of the 
dead-time counter and is explained in greater detail in 
Appendix A. 

Note that the dead time, during which the counter 
is unable to record any incoming photons, induces sat- 
uration of the maximally possible counts to JWmax ~ 
limp_j,tx) C(p) = 1/S. As mentioned before, the stan- 
dard gluing procedure involves, before fitting the ana- 
log and photon counting traces, correcting the counts m 
for the dead-time effects with the inverse of the function 
in Eq. (3), 



p = C ^{m) = 



m 



\ — 5m' 



(5) 



Unfortunately, the inverse has a singularity aim = JWmax 
and produces negative photon estimates p for m > 
ntmax- In the standard gluing procedure this, and the 
fact that the dead time t is not well known, limits the 
range of useful data of the measured photon counting 
traces. 

As we will see later on, the procedure developed here 
does not suffer from this drawbacks since only the non- 
singular function C(p) is used and the estimation of the 
dead-time value is part of the method. 



C. Overlap region 

In the case of a large number of incoming photons, 
p 2> 1/5, only the analog signal carries useful infor- 
mation due to the inevitable saturation of the photon 
counter. For a small photon influx the situation is re- 
versed since the analog signal has reached the levels of 
the electronic noise w^hile the photon-counting is in the 



^ The procedure given here can be naturally adapted also for the ex- 
tending (or paralyzable) type of photon counters by replacing Eq. (3) 
with C(p) = pexp(i5p) and its associated variance [2], 



ideal proportional mode with almost no dead-time ef- 
fects. Therefore, outside of the relative overlap region, 
the quality of data of one or the other mode prevails. 
From the simple measurement models given above we 
thus require good accuracy in the overlap region and 
that the winning model is supplying the correct solution 
far away from the overlap. 

Due to these relaxed requirements we can approx- 
imate the variance of the dead-time affected photon 
counter in Eq. (4) with the Poisson variance of the 
photon-count itself (see Appendix A), 



V\m] 



m 



dp). 



(6) 



D. Summation of lidar traces 



It is quite common practice in lidar measurements to 
additionally increase the dynamical range of the data ac- 
quisition by summation {time stacking) of consecutive li- 
dar returns. With fast laser-pulse repetition rates it is 
reasonable to assume that within the summation time 
the atmosphere does not introduce substantial sources 
of additional variance beyond the natural Poisson-like 
fluctuations of the backscattered photons. 

Denoting by fls the sum of Ng analog measurements a 
at the same range of consecutive lidar traces and with p^ 
the sum of the arrived photons, the photon conversion 
in Eq. (1) is transformed into flg = NsA{ps/Ns) = aps + 
Ns/5. The variance in Eq. (2) scales as V[fls] = ^sl^- 

The mean photon count obtained by summation of Ng 
consecutive lidar returns, NsC{ps/Ns), has a nice prop- 
erty of retaining the general form of Eq. (3) with the 
dead-time fraction effectively transformed into 5/Ns- 
Nevertheless, for large photon numbers p the variance 
depends on summation in a non-linear way and has to 
be evaluated as NsV^(ps/Ns). 

The measurement models given above are thus, at 
least to some degree, invariant with respect to the sum- 
mation as long as the following transformation of the 
acquisition parameters is taken into account. 



a^ a, ^ H^ ^/Ns, 7^ ^ J^/N^, 



5 ^ Ns5. (7) 



III. INITIAL ESTIMATES 

For any nonlinear minimization procedure it is of ut- 
most importance to acquire accurate initial values for 
minimized parameters. In our case the initial values 
for parameters a, j6, 7^ are obtained from a least-squares 
minimization of 



Amin 






(8) 



where only analog and photon-counting data points 
{ai,mj) of the lower left corner are used (see Fig. 1), i.e. 
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Figure 1. Plot of analog vs. photon-counting data points 
{aj,mj) for Ng = 20 summed lidar returns. The slanted line 
on the left is a fit for the photon-to-analog conversion parame- 
ters (X, p and 7^ in the range of the lower left corner (indicated 
by the left arrow-box) with dashed lines illustrating variance 
a ± 2j. The horizontal line is a fit for dead-time fraction S in the 
range of the upper right corner (right arrow-box). The thick 
line is the resulting prediction for m = C{A^^ (a)). 



the lower 10% of the whole photon-count range. Ini- 
tial estimates for the "gluing" parameters a. and /5 are 
thus obtained in a manner similar to the standard pro- 
cedure suggested by the manufacturers of the transient 
recorders [6] or other more detailed studies [5]. 

Since the fluctuations of the analog data in the lower 
left corner of Fig. 1 are dominated by the electronic 
noise, an estimate for the analog variance 7^ is obtained 
simply from the residuals in Eq. (8), 7^ = Xmin/^dof' 
where the number of degrees of freedom is N^gf = N — 2 
with N the number of data points involved in the fit. 
From this point on, 7^^ is kept fixed and will not be the 
subject of minimization. 

Fitting the photon counts m to a constant in the tail 
of the large analog values (upper right corner of Fig. 1) 
gives an estimate for the dead-time fraction, S = 1/ {m). 
For the {m) estimate typically only the data in the largest 
30% of the analog values has been used, but excluding 
all data points with ADC saturation (which we discard 
in this procedure anjrway). 

See Fig. 1 for the results of the initial fits to an example 
lidar return which will be used throughout this analysis. 
The data was obtained with a back-scatter lidar at wave- 
length of 355 nm, pulse repetition rate of 20 Hz and trace 
summation Ns = 20. The light sensor was a Hamma- 
matsu R3200 photomultiplier tube connected to a high- 
voltage of approximately 800 V. The return signal was 
acquired with a Licel TR 40-160 transient recorder oper- 
ated at a sampling frequency of 40 MHz. The recorder 
is delivering discretized analog signal traces with 12 bit 
ADC resolution and photon-counting traces with a max- 
imal count rate of 250 MHz, both with trace depth of 16k 



samples [7]. The example trace was recorded with 20° 
elevation on a relatively clear night and contains only a 
thin and faint layer of haze around the range of 13.5 km. 



IV. MAXIMUM LIKELIHOOD 

From the two measurement models described above 
we can construct a likelihood C for the total trace as a 
product over all trace time elements i of a likelihood of 
observing p, photons given the analog measurement fl, 
and the photon count rrii, 



C =Y[C{ai,mi,pi), 



(9) 



where likelihood C{aj,mj, pj) is a product of the proba- 
bility P(fl/|p/) to observe an analog signal and the prob- 
ability to have a certain photon count P(ot,|p,) given a 
number of photons. We can model the analog proba- 
bility with the normal (Gauss) distribution J\f{x,cr^) = 
ex-p{—x^/2a^)/V2TW^ using the linear transformation 
Eq. (1) and the corresponding variance 7^. According 
to Eq. (6), the photon count probability can be approx- 
imated with the Poisson distribution from Eq. (Al) so 
that the resulting likelihood is expressed as 

C{ai,mi,pi) = M{ai - A{pi),j^) x P„,,(C(p,)) (10) 

The corresponding deviance is defined as 

D = -21n£ = -2j2^nC{ai,mi,pi) = 

i 

= J2V{ai,mi,pi), 



(11) 



where 

T>{ai,mi,pi) = In27r7^ 



72 



+ 2[\nmi\ + C{pi) - m,\nC{pi)] (12) 

is the deviance for a particular data point . The motiva- 
tion for using the deviance version of likelihood comes 
from the fact that for the normal-like distribution prob- 
abilities the deviance is equivalent to the usual x^ es- 
timator. Nevertheless, the minimization of Poisson-like 
distribution probabilities can not be formulated in terms 
of a simple x^ formalism. 

The solution to the minimal deviance (or, in other 
words, maximal likelihood) 



T) = min. 
is usually found by solving for an extreme 

VD = 0, 

^ with specific requirement that x In = 



(13) 
(14) 



where the gradient V is formed by derivatives over 
the whole parameter space. In our two-measurement 
model, the deviance (likelihood) depends on the follow- 
ing parameters: the a. and /5 coefficients from the analog- 
to-digital conversion A{p), the variance of the analog 
signal 7^, and the dead-time fraction 3 of the photon 
counter. In addition to these four model parameters, the 
deviance depends also on all (unknown) numbers of in- 
coming photons Pi- In general, the deviance thus has 
N + 4 parameters for 2N data points (analog and pho- 
ton counts). 

Due to the particular structure of the deviance in 
Eq. (11), we can split the global minimization procedure 
for I) in two parts: the outer part drives the minimiza- 
tion over a, fi, j'^ and S parameters, and the inner part 
deals with the "nuisance" parameters p, for each itera- 
tion of the outer part. 

Assuming that the outer part already supplies pa- 
rameters a, j6, 7^ and 3, the inner part proceeds as fol- 
lows: since in the total deviance V only the fth term 
'D{ai,mi,pi) depends on particular p,, we can simplify 
the Pi part of its gradient. 



dV dT>[ai,mi,pi) 
dpi dpi 



(15) 



by introducing a marginalized (conditional) number of 
photons. 



pi = ar^mh\'D{ai,mi,pi), 
Pi 

and profile deviance (as in profile likelihood) 



(16) 



'D{ai,mi) = min'D(ai,mi,pi) = T>(ai,mi,pi). (17) 

Pi 

Solving this equation for all p, produces a total deviance 
without the nuisance parameters p,. Finally, the de- 
viance is contracted into 



V = Y^t>{ai,mi), 



(18) 



w^hich depends only on the remaining four parameters 
^, jS/ 7^/ and 3, and is minimized by the outer part of the 
procedure. 

For our two particular measurement models, 
Eqs. (15)-(17) in the inner part of the minimization 
correspond to finding a suitable root of a polynomial 
of the fourth order in p,-. Although analytical solutions 
exist, they are not very practical for real application. 
The solution can be obtained with the Newton-Raphson 
iteration 



p[;-+i] = p\i] _ '^'i.ai,mi,pf) 
' V"{a„mi,pf) 



(19) 



where V and V" are respectively the first and the 
second derivative of 'D{ai,mi,pi) with respect to p. 
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Figure 2. Dependence of the normalized deviance V/N on 
relative offset toffset between the analog and photon-counting 
traces. The arrow indicates the position of the minimum at 
^offset =4A< = 100 ns. 



in Eq. (12). The iteration in Eq. (19) is started with 

a suitable approximation pj ' and is repeated until 

I ~ Pi becomes smaller than e, with e set to some 

small number {^ 10^^). Note that in some cases the 
minimum over p, is not the zero of the derivative in 
Eq. (15), but can be instead found at the boundary, p, = 
0, of the validity interval of the p, parameter 

The outer part of the minimization deals with the total 
deviance in Eq. (18) with respect to the remaining non- 
fixed parameters'' and this can be carried out with a va- 
riety of nonlinear minimization procedures (see for ex- 
ample [8]). Denoting the final values of the parameters 
in the deviance minimum with a, ^, and 3, the set of final 
values of nuisance parameters 



Pi{^,^,l^J) 



(20) 



in the global minimum of the deviance represents the ul- 
timate (most likely) synthesis of the analog and photon- 
counting modes of the lidar data acquisition. Note that 
7^ is kept fixed at the value of the initial approximation 
throughout this procedure. 



A. Relative acquisition delay 

In the acquisition system the input signal flows 
through quite different electronic sub-components of 
the transient recorders (e.g. see schematic in manual 
[7]) so expecting hardware and firmware related differ- 
ences in delay time is highly justified. The final ana- 
log and photon counting traces can thus be subject to 



^ The minimization in Eq. (17) is thus embedded inside the outer min- 
imization. 



70 


1 1 1 1 1 1 1 


i 


60 




fiih 


50 






40 






30 






20 




1 ' 


10 


:W' 










-— analog 

^ shifted analog 

^ photon counts 




^/^'"^; 



3540 3560 3580 3600 3620 3640 3660 
I 

Figure 3. Haze feature before and after the shift of the 
analog trace. The interval [3450,3660] in i corresponds to 
[13.275, 13.725] km in range. 
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substantial relative offset. In the case of our particu- 
lar recorder this offset amounts to four sample times, 
^offset = 4 At = 100 ns. The dependence of the total de- 
viance on ^offset is shown in Fig. 2. How this offset is 
influencing the details of the example trace can be ob- 
served in Fig. 3. The two haze features found around 
13.5 km in the analog and (uncorrected) photon-count 
modes perfectly match after the foffset shift. The same 
holds for the small noise-like features in the rest of the 
trace, mostly responsible for the distinct minimum of 
deviance in Fig. 2. 



B. Parameter bias due to data distribution 

Typical lidar returns do not cover uniformly the 
whole dynamic range available in analog, a, and 
photon-counting, m, modes. Most of the data resides 
in the tail of the lidar return, occupying only the lower- 
left sector of the (a, m) phase space (c.f. Fig. 1). The 
large contribution of this sector can thus introduce a 
bias in the likelihood maximization procedure. Further- 
more, the acquisition parameters are influenced by the 
different parts of the {a,m) phase space. Analog base- 
line /5 is well defined by the lidar tail (lower left sec- 
tor in Fig. 1), the dead-time fraction 6 is mostly sensi- 
tive to the large-signal parts (upper right sector), and 
the photon-to-analog coefficient a is mostly influenced 
by the small and medium part of the trace (left side in 
Fig. 1). To remove and quantify this bias we can bin 
the data with several different partitions of the {a,m) 
phase space and balance the relative contribution of 
each data point {ai,mj) to the total deviance with a 
weight w{ai,mi) — Wj, where / is the appropriate bin 
index. To maintain the correspondence with the pre- 
vious non-binned case (equivalent to the binning with 
Wj = 1) all binning variants are required to satisfy a 
summation rule X]/ Wj = N, with N being the number 



Figure 4. Example of a debiasing binning of data points {aj,mj) 
from the example lidar return into the fan-like bins. Note that 
the bottom bin contains many more data points than all of the 
other bins. The total deviance weights are set proportionally to 
the inverse of the particular bin count N,, which is also shown 
for all bins. 



of data points. In all binning variants the weight should 
be proportional to the inverse density so that all points 
in a particular bin contribute to the deviance with the 
same weight as the other non-empty bins. With these 
two requirements the weights are obtained as 



w{ai,mi) = Wj 



N 



Ne>Nj' 



(21) 



where N^ is the number of non-empty bins and N, the 
number of data points in a particular bin j. 

The debiased version of deviance from Eq. (18) is now 
written as 



V = Y^w{ai,mi)V{ai,mi). 



(22) 



Here we consider several variants of the binning divi- 
sions. 

• Un-binned case, w(ai,mi) = 1: every data point 
is considered with the same weight. Note that in 
this way in cases of lidar traces with long tails, 
the overwhelming contribution will come from the 
points with small values in both modes fl, and m,. 

• Fan-like binning: data is categorized into a his- 
togram with fan-like bin shapes radiating from the 
lower right corner, (flmax/Wmin) = (2^2,0) (see 
Fig. 4-left). 

• Infinitely-fine binning: since both measurement 
modes a, and m, have only discrete integer val- 
ues, we treat each of these discrete points with the 
same weight. 
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Figure 5. Dependence of the analog parameters a, ji and 
photon-count parameter S on the number of non-empty bins 
Ng for the fan-like binning (black points). Left and right ar- 
rows indicate parameter values from the un-binned and the 
infinitely-fine binned cases, respectively. The values from the 
initial estimates are shown in dashed and dotted lines (see the 
text for more details). 



In the top two panels of Fig. 5 initial estimates for 
the analog parameters ec and /5 are shown for two ver- 
sions of the x^ fit from Eq. (8). The dashed lines are 
for the x^ with squares of a, — A{mi) (involving uncor- 
rected photon-counts) and the dotted lines are for a ver- 
sion where the dead-time fraction S is estimated first and 
then the x^ is constructed with dead-time corrected pho- 
ton counting, a, — A(C~^(wi,)). In the bottom panel of 
Fig. 5 the initial estimate for dead time S is shown as a 
dashed line. The arrows on the left and the right of all 
panels are indicating the results of the un-binned and 
the infinitely-fine binned cases, respectively. The rest of 
the points correspond to the fan-like binning for differ- 
ent granularity of the binning represented by the num- 



ber of non-empty bins Nq- 

While there are only minor differences between the 
un-binned and infinitely binned results, the fan-like bin- 
ning exhibits large variations of the order of 25% for 
a, 5% for /5, and 15% for S w^hen changing the binning 
size"*. Nevertheless, with increasingly fine binning the 
parameters tend to converge to the two values obtained 
by the un-binned and infinitely binned cases, indicat- 
ing that (at least for this trace depth) the maximum like- 
lihood method is only mildly biased by the particular 
data distribution. On the other hand, in longer traces 
the distribution of the data in the lower left sector might 
influence the reconstruction, especially if the measure- 
ment models are not accurate enough. 



V. RESULTS 

A. Reconstructed number of photons 

In order to follow the evolution of the reconstructed 
(most likely) number of photons p from Eq. (20) let us 
introduce a transition indicator 



(23) 



where pa and pm are direct estimates for photons ob- 
tained from the two measurements. 



Oa=A-\a) 



a- ji 



(24) 



is the analog-to-photon conversion (inverse of Eq. (1)) 
and 



Pm = C ^{m) = - 



m 



5m 



(25) 



is the corrected photon count (inverse of Eq. (3)). Note 
that the latter works only for m < 1/S. The values of 
the indicator u defined in this way will be close to 1 
when the reconstructed number of photons p is close 
to the prediction from the current inversion pa- On the 
other hand, u will be close to when p is close to the 
dead-time corrected photon-count p,„. Values between 
and 1 are indicating transition between the two mea- 
surements. 

In Fig. 6 the changing of the value of indicator u along 
our example trace is shown. We can clearly identify 
three regions of the indicator 's behavior. 

For indices i below ~ 2000 the indicator m w 1 and 
thus the reconstructed number of photons is closely fol- 
lowing the estimate from the analog signal. In this re- 
gion the photon counting is saturated and, due to the 



* Uncertainties are in fact not so large, considering that the parame- 
ters are obtained on a single trace with Nj = 20 summation only. 




Figure 6. Behavior of the transition indicator u from 
Eq. (23). The [1500,5000] interval in i corresponds to the 
[5.625, 18.75] km interval in range. 



dead-time, its resolution is heavily suppressed. The 
maximum likelihood method thus shifts the result to- 
wards the more accurate measurement: the analog sig- 
nal. 

For indices i above ^ 4000 the indicator m ~ and the 
reconstructed number of photons is closely following 
the photon-counting. Here the analog signal is diving 
into the noisy region around the analog baseline where 
SNR is small. On the other hand, the photon counting is 
far from being effected by the dead-time and thus max- 
imum likelihood gives it deserved emphasis. 

In the intermediate region, for indices i between -^ 
2000 and ^ 4000 where < m < 1, the reconstructed 
number of photons lies between the two measurements 
w^hich both have degraded accuracy, analog signal due 
to poor SNR and photon-counting due to the dead-time 
effects. Nevertheless, the value of p is chosen according 
to the maximum likelihood, effectively combining the 
two less accurate measurements into one with smaller 
error^. 

Note that all variants of the standard "gluing" pro- 
cedures would in this picture produce a step-like func- 
tional form of indicator u, abruptly crossing over from 1 
to at a position that depends on a particular choice of 
the "gluing" method. 



B. Multiple lidar traces and dynamic range 

Standard processing of the lidar returns usually em- 
ploys extensive stacking (summation) of the lidar traces 
in order to increase the SNR ratio. Since the maximum 
likelihood method has no built-in notion of range and 



data-point ordering, multiple traces can instead be con- 
catenated in order to increase the accuracy of the recon- 
struction. All points from different traces can thus be 
treated equally and processed at the same time, as long 
as the acquisition parameters a, j5, 7^, and 3 are consid- 
ered stable in the respective time frame of the recording 
of the traces. Nevertheless, this procedure will suffer a 
slowdown linear with the number of data points used'^. 

In case of our equipment, relative scatter of the recon- 
structed parameters within individual 480 s runs were 
below 1.6% for a, 0.24% for (S, and 0.28% for S. Rel- 
ative differences between the mean values of the re- 
constructed parameters for the runs from the beginning 
and the end of the measurement campaign were below 
0.61% for a, 0.54% for ^, 1.62% for 7^, and 0.18% for S. 
In such stable conditions it is thus possible to concate- 
nate a large number of lidar traces in order to increase 
the accuracy of the reconstructed photon returns. 

In the usual lidar operation, the gain of the analog 
channel is set to a level which produces a discretized 
trace with the maximum signal slightly below (or, like 
in our case, above) the ADC saturation limit. In this way 
the analog signal can cover the whole dynamic range of 
the ADC, exposing nicely the saturation of the photon 
counting (see Fig. 1). It turns out that if for whatever 
reason the analog signal is covering only a smaller frac- 
tion of the available dynamic range or that the amount 
of photons is not saturating the photon counting, the 
procedure described above will still produce stable and 
reasonable estimates of the a, ji, and 7^ parameters since 
they are anyway dominated by the data in the lower- 
left corner of the («;,;«,) diagram. On the other hand, 
the dead-time fraction S will in this case tend to zero 
(towards the ideal counter), since the absence of the 
photon-counting saturation in the data effectively gives 
an estimate 1 / .5 — )• 00, but all this without actually influ- 
encing the reconstructed number of photons p,. 



VI. DISCUSSION AND CONCLUSIONS 

Modem transient recorders offer two principally dif- 
ferent measurements of the same photon return: the dig- 
itization of the analog signal and the photon-counting 
mode. We are thus challenged, not only to use them in 
their respective regions of validity (like the usual "glu- 
ing" methods) but to combine them into a more accurate 
estimate of the photon numbers by using detailed statis- 
tical models of both acquisition processes. 

The maximum likelihood procedure described in this 
work offers a reconstruction of photon returns that has 
a natural transition between the analog and photon- 
counting signals and is based on their analytical mea- 
surement models. In this work we have been using 



^ e.g. maximum-likelihood combination of two normally distributed 
measurements with errors Uj and 02 gives a new estimate with a 



smaller error criO'2 1 y (j\ + ov 



^ At the time of the writing, it takes 0.2 s per 16k trace on a normal 
desktop computer. 
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fairly rudimentary models of the two measurement pro- 
cesses, nevertheless they still adequately capture the 
main strengths of the new reconstruction procedure. 
Furthermore, if more detailed models are needed, they 
can be simply included into the probability expressions 
entering the likelihood function. 

In this method we are strongly discouraged to attempt 
any kind of background subtraction. The background 
photons are treated as a normal signal since they ap- 
pear in both measurement modes as viable data. Any 
removal of background (from dawn or daylight) should 
be done on the final reconstructed photon numbers. 

The maximum likelihood method works in all condi- 
tions, even in the presence of clouds or other enhanced 
scattering objects. It fails only when the input levels of 
photons are not exploring the whole dynamic range of 
the transient recorder (i.e. both comers of Fig. 1) and 
thus no reliable estimate of the acquisition parameters 
a, /3, 7^, and 5 can be obtained. Through the offset anal- 
ysis of the likelihood value it offers a simple way for es- 
timation of the potential relative delay between the two 
measurement traces, which in our particular case turned 
out not to be negligible. 

The code implementing the maximum-likelihood re- 
construction of lidar returns is available under the GPL3 
license at http://www.ung.si/~darko/lidar/. 
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Appendix A: Dead-time counter 

For an ideal counter with sampling time Af the dis- 
crete probability distribution of the number of counts k 
for a Poisson process with rate r and mean M = rAf is 
given by "Pj-C-^) "where 



V^{x) 



X e ^ 
k\ ■ 



(Al) 



For a counter with non-extending dead-time t [9] the 
corresponding probability distribution can be found in 
Ref. [10]. Restructuring the equations and performing 
partial summations of the expressions given there, the 
probability Wj. of observing k counts now becomes 



Wk 



1 



1+M3 



[K/t-i - 2Rk + Rk+i + Afc] , 



(A2) 



where S = r/Af, the short-hand R^ = Rkih) ^nd the 
truncated mean for k counts is given by 



tk = M{l-kS). 
^/c = J^kitk) is fully expressed as 



(A3) 



k-l 

Rk{x) = U{x)J^{k-j)Vj{x) = 
j=o 

= U{x)[{k-x)Q{k,x)+kVk{x)], (A4) 

where Q{i,x) = T{i,x)/T{i) is the regularized upper in- 
complete Gamma function [11], with the upper incom- 
plete Gamma function 



T{a,x) = I u''^'^e~"du (A5) 



and r(fl) = r(fl, 0). The unitary step function U{x) is 
defined in the usual way. 



U(x) 



1 if X > 0, 
otherwise. 



(A6) 



and the remainder A^ in Eq. (A2) is given explicitly as 

'O iffc^X-l, 

Ak= {{K + 1){1+MS)-M iik^K, (A7) 

M-K{1+MS) iik = K + l. 

The upper limit on possible counts depends on the 
dead-time fraction. 



K= [1/S\, 



(A8) 



where [x\ is denoting the largest integer smaller (and 
not equal) than x. 

The mean dead-time count m can be expressed in 
terms of the ideal count M as 



m = C{M) 



M 
1+MS' 



(A9) 



The exact expression for the variance of the dead-time 
counter is 



K 



V, 



+ H{m-K), 



,Ei(^-ik)Qik,tk)+m{tk)] + 



(AlO) 



where the "hump" function is defined as H{x) = x{l — 
x). 

Using d = MS as a mean "lost" count, for d ^ 1 
the exact variance from Eq. (AlO) asymptotically [10] be- 
haves as 



y,«M 



1 ^^{6+4:d + d^) 



(l+d)3 6Mil+dy 



(All) 



where the part in the square brackets is the suppression 
factor relative to the Poisson variance Vpoisson = A/f . For 
d 2> 1 the variance is well described by a fully saturated 
dead-time counter, 

Vs w H{m -K)= H{hac{m)), (A12) 

where frac(x) = x — [xj is the function returning frac- 



tional (non-integer) part of an argument. Note that the 
expression for asymptotic variance in Eq. (All) con- 
verges in this regime to 1/6 and thus well describes 
the average of the oscillatory dependence of V^ on S in 
Eqs. (AlO) and (A12). 
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